//
//  KCurve.m
//  Kia
//
//  Created by Andrey on 03/04/2009.
//  Copyright 2009 Karma Software. All rights reserved.
//

#import "KCurve.h"

@implementation KCurve

NSInteger ComparePoints(id p1, id p2, void* context)
{
	NSPoint point1 = [p1 pointValue];
	NSPoint point2 = [p2 pointValue];
	
	if (point1.x < point2.x)
		return NSOrderedAscending;
	else if (point1.x > point2.x)
		return NSOrderedDescending;
	else
		return NSOrderedSame;
}

@synthesize controlPoints;

- (id) init
{
	if (self = [super init])
	{
		controlPoints = [[NSMutableArray arrayWithCapacity:MAX_CURVE_CONTROL_POINT_COUNT] retain];
		splines = (CubicSpline*)malloc(sizeof(CubicSpline) * (MAX_CURVE_CONTROL_POINT_COUNT - 1));
		needRecomputeCurve = YES;
	}
	
	return self;
}

- (void) addControlPoint: (NSPoint) point
{
	[controlPoints addObject:[NSValue valueWithPoint:point]];
	[controlPoints sortUsingFunction:ComparePoints context:nil];
	needRecomputeCurve = YES;
}

- (void) replaceControlPointAt: (NSUInteger)index withPoint: (NSPoint)point
{
	[controlPoints replaceObjectAtIndex:index withObject:[NSValue valueWithPoint:point]];
	[controlPoints sortUsingFunction:ComparePoints context:nil];
	needRecomputeCurve = YES;
}

- (void) removeControlPointAt: (NSUInteger)index
{
	[controlPoints removeObjectAtIndex:index];
	needRecomputeCurve = YES;
}

- (void) removeControlPoint: (NSPoint) point
{
	[controlPoints removeObject:[NSValue valueWithPoint:point]];
	needRecomputeCurve = YES;
}

- (void) removeAllControlPoints
{
	[controlPoints removeAllObjects];
	needRecomputeCurve = YES;
}

- (void) recomputeCurve
{
	if (needRecomputeCurve == YES)
	{
		for (int i = 0; i < [controlPoints count]; i++)
		{				
			splines[i].a = [[controlPoints objectAtIndex:i] pointValue].y;
		}
		
		CGFloat* steps = 
		(CGFloat*)malloc(sizeof(CGFloat) * ([controlPoints count] - 1));
		for (int i = 0; i < [controlPoints count] - 1; i++) 
		{
			steps[i] = 
			[[controlPoints objectAtIndex:(i + 1)] pointValue].x - 
			[[controlPoints objectAtIndex:i] pointValue].x;
		}
		
		CGFloat* alpha = 
		(CGFloat*)malloc(sizeof(CGFloat) * ([controlPoints count] - 1));
		for (int i = 1; i < [controlPoints count] - 1; i++) 
		{
			alpha[i] = 
			3 * (splines[i + 1].a - splines[i].a) / steps[i] -
			3 * (splines[i].a - splines[i - 1].a) / steps[i - 1];
		}
		
		CGFloat* l = 
		(CGFloat*)malloc(sizeof(CGFloat) * [controlPoints count]);
		CGFloat* m = 
		(CGFloat*)malloc(sizeof(CGFloat) * [controlPoints count]);
		CGFloat* z = 
		(CGFloat*)malloc(sizeof(CGFloat) * [controlPoints count]);
		l[0] = 1.0;
		m[0] = z[0] = 0.0;
		for (int i = 1; i < [controlPoints count] - 1; i++) 
		{
			l[i] = 
			2 * ([[controlPoints objectAtIndex:(i + 1)] pointValue].x - 
				[[controlPoints objectAtIndex:(i - 1)] pointValue].x) - 
				steps[i - 1] * m[i - 1];
			m[i] = steps[i] / l[i];
			z[i] = (alpha[i] - steps[i - 1] * z[i - 1]) / l[i];
		}
		l[[controlPoints count] - 1] = 1.0;
		z[[controlPoints count] - 1] = 0;
		splines[[controlPoints count] - 1].c = 0.0;
		
		for (int i = [controlPoints count] - 2; i >= 0; i--)
		{
			splines[i].c = z[i] - m[i] * splines[i + 1].c;
			splines[i].b = 
				(splines[i + 1].a - splines[i].a) / steps[i] - 
				steps[i] * (splines[i + 1].c + 2 * splines[i].c) / 3;
			splines[i].d = (splines[i + 1].c - splines[i].c) / (3 * steps[i]);
			splines[i].t0 = [[controlPoints objectAtIndex:i] pointValue].x;
/*			
			NSLog(@"Spline #%d: {%6.4f, %6.4f, %6.4f, %6.4f, %6.4f}", 
				  i, splines[i].a, splines[i].b, splines[i].c, splines[i].d, splines[i].t0);
*/
		}
		
		// Cleanup
		free(z);
		free(m);
		free(l);
		free(alpha);
		free(steps);
		
		needRecomputeCurve = NO;
	}
}

- (CGFloat) valueAt: (CGFloat)t
{	
	if (t <= 0)
		return 0.0;
	
	if (t >= 1.0)
		return 1.0;
	
	int interval = 0;
	
	while (interval < [controlPoints count] && t > [[controlPoints objectAtIndex:interval] pointValue].x) 
	{
		interval++;
	}

	if (needRecomputeCurve == YES)
		[self recomputeCurve];
	
	CubicSpline s = splines[interval - 1];
	t -= s.t0;
	
	CGFloat value = s.a + s.b * t + s.c * t * t + s.d * t * t * t;
	
	if (value < 0.0)
		value = 0.0;
	else if (value > 1.0)
		value = 1.0;
	
	return value;
}

#define CURVES_WINDOW_SIZE 260
- (void) saveToTikzFile: (NSString*)filename
{
	NSString* imageName = 
	[filename stringByDeletingPathExtension];
	
	NSString* plotDataFilename = 
	[imageName stringByAppendingString:@".table"];
	
	NSString* plotData = [NSString string];
	for (int i = 0; i < CURVES_WINDOW_SIZE; i++)
	{
		CGFloat x = 2 * (CGFloat)i / CURVES_WINDOW_SIZE;
		CGFloat y = 2 * [self valueAt:(0.5 * x)];
		
		plotData = 
		[plotData stringByAppendingFormat:@"%.5f\t%.5f\n", x, y];
	}
	[plotData writeToFile:plotDataFilename 
			   atomically:NO 
				 encoding:NSUTF8StringEncoding 
					error:NULL];
	
	NSString* tikzContent = @"\\begin{tikzpicture}[scale = 3, >=stealth]\n";
	// Draw grid
	tikzContent = 
	[tikzContent stringByAppendingString:@"\t\\draw[very thin, color = gray] (-0.1, -0.1) grid (2.1, 2.1);\n"];
	// Draw axis
	tikzContent = 
	[tikzContent stringByAppendingString:@"\t\\draw[->] (-0.1, 0) -- (2.1, 0) node[right] {$c_b$};\n"];
	tikzContent = 
	[tikzContent stringByAppendingString:@"\t\\draw[->] (0, -0.1) -- (0, 2.1) node[above] {$c_b^*$};\n"];
	tikzContent = 
	[tikzContent stringByAppendingFormat:@"\t\\draw[thick] plot[smooth] file {%@};\n", plotDataFilename];
	tikzContent = 
	[tikzContent stringByAppendingFormat:@"\t\\draw (0, 0) -- (2, 2);\n"];
	
	// Draw control points
	for (int i = 1; i < [controlPoints count] - 1; i++)
	{
		NSPoint controlPoint = [[controlPoints objectAtIndex:i] pointValue];
		
		tikzContent = 
		[tikzContent stringByAppendingFormat:@"\t\\draw (%.4f, %.4f) circle (0.01cm);\n", 
		 2 * controlPoint.x, 2 * controlPoint.y];
		
		tikzContent = 
		[tikzContent stringByAppendingFormat:@"\t\\draw[<->] (%.4f, %.4f) -- (%.4f, %.4f);\n", 
		 2 * controlPoint.x, 2 * controlPoint.x, 
		 2 * controlPoint.x, 2 * controlPoint.y];
	}
	
	tikzContent = [tikzContent stringByAppendingString:@"\\end{tikzpicture}"];
	[tikzContent writeToFile:filename 
				  atomically:NO 
					encoding:NSUTF8StringEncoding 
					   error:NULL];
}

- (void) dealloc
{
	free(splines);
	[controlPoints release];
	[super dealloc];
}
@end
